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ABSTRACT 

We present the results of a new, non-parametric method to reconstruct the 
Galactic dark matter profile directly from observations. Using the latest kine¬ 
matic data to track the total gravitational potential and the observed distribution 
of stars and gas to set the baryonic component, we infer the dark matter con¬ 
tribution to the circular velocity across the Galaxy. The radial derivative of this 
dynamical contribution is then estimated to extract the dark matter profile. The 
innovative feature of our approach is that it makes no assumption on the func¬ 
tional form nor shape of the profile, thus allowing for a clean determination with 
no theoretical bias. We illustrate the power of the method by constraining the 
spherical dark matter profile between 2.5 and 25kpc away from the Galactic cen¬ 
tre. The results show that the proposed method, free of widely used assumptions, 
can already be applied to pinpoint the dark matter distribution in the Milky Way 
with competitive accuracy, and paves the way for future developments. 
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1. Introduction 


Mapping out the distribution of dark matter in our Galaxy rests as a paramount task 
with potentially far reaching implications for astroparticle physics and cosmology. This 
is important to understand galaxy formation and to feed searches aimed at unveiling the 
very nature of dark matter. In particular, direct and indirect dark matter searches rely 
heavily on the findings of numerical simulations. It is therefore essential to extract the 
Galactic dark matter distribution directly from observations. In the outer Milky Way (at 
Galactocentric radii greater than ~ 20kpc), where baryons contribute little to the total 
mass budget, the gravitational potential traces closely the dark matter component and the 


total mass enclosed can be constrained using convenient tracers (e.g., Sakamoto et al. 2003 


Dehnen et al. |2006 Xue et al. 2008 Bhattacharjee et al. 2014 |Kafle et al. 2014), although 


with important degeneracies in the tracer population modelling. By contrast, in the inner 
Galaxy (i.e., in the inner ~ 20kpc) the baryonic contribution is very significant and its 
morphology rather uncertain, which makes the evidence for dark matter difficult to establish 


and the extraction of its distribution a delicate undertaking ( 

Iocco et al. 

2015 

). This has 

been addressed by many authors with different methods (e.g., 

Dehnen & Binney|l998; 

Sofue 

et al. 2009 Catena & LUlio 2010; Weber & de Boer 2010; Iocco et al. 2011 

Nesti & Salucci 

2013 

Bovy & Rix 2013 

Loebman et al. 

2014 

), all of which do, however, make explicit 


assumptions about the underlying dark matter profile: typically, a multi-parameter profile 
is fitted to the observations together with a given baryonic component. The class of “local” 
methods to measure the dark matter density in the solar neighbourhood (e.g., Salucci et al. 


2010 Smith et al. 2012 Garbari et al. 2012 Zhang et al. 2013; Read 2014) avoids this bias, 


yet such methods are not easily applicable elsewhere in the Galaxy. An approach free of 


profile assumptions has been developed and successfully tested in external galaxies (Persic 


et al. 1996; Salucci et al. 2007), but never applied to our own Galaxy given the sizeable 


uncertainties of both kinematic data and baryonic modelling. 


In this Letter we show that the latest rotation curve measurements and baryonic mod¬ 
els make it possible to infer the dark matter profile directly from Milky Way observations 
without unnecessary assumptions. Our results confirm that this approach is quantitatively 
competitive to the others used so far, while presenting the noticeable advantage of making 
no a priori assumption on how dark matter is distributed across the Galaxy. 


2. Methodology 

The total gravitational potential of our Galaxy can be written as a sum of two compo¬ 
nents, namely baryons and dark matter: (j) to t = (ftb + 4>dm- The left-hand side (or rather its 
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radial derivative) is fixed by measurements of the rotation curve, whereas the first term on 
the right-hand side is set by the observed distribution of stars and gas. These are the two 
data inputs needed to infer the distribution of Galactic dark matter. 

Regarding the rotation curve, we determine the angular circular velocity u> c with a broad 
collection of tracers comprising gas kinematics (HI and CO terminal velocities, HI thickness, 


HII regions, giant molecular clouds; Fich et al. (1989); McClure-Griffiths & Dickey (2007); 


Luna et al. (2006); Honma & Sofue 

(1997) 

; Brand & Blitz ( 

1993); Hou et al. 

2009)), star 

kinematics (open clusters, planetary nebulae, classical cepheids, carbon stars; 

Frinchaboy 

& Majewski (2008p; Durand et al. i 

1998); 

Pont et al. (1997 

); Battinelli et al. 

(2013)) and 


masers (Reid et al. 2014) in a total of 2780 measurements distributed across Galactocentric 


radii R = 0.5 — 25kpc. Our compilation of data improves upon commonly used compilations 


(Sofue et al. 2009 Bhattacharjee et al. 2014) by including numerous tracers available in 


the literature but often neglected. For each object, the measured line-of-sight velocity in 
the local standard of rest v\™ is converted into the angular circular velocity through v\°J, = 
(RqUJc ~ Vo) cos b sin £, where £, b are the Galactic longitude and latitude, R 0 is the distance 
to the Galactic centre and u 0 = v c (Rq) the local circular velocity. Throughout the analysis 


we take Rq = 8kpc, Vq = 230km/s and the peculiar solar motion of Schonrich et al. (2010). 


The uncertainties on both distance and kinematics are assigned according to each source 
reference and propagated to R and oj c , respectively. We have checked that our determination 


of the rotation curve is solid against systematics due to spiral arms (Brand & Blitz 1993) 


and against the non-circularity of tracer orbits. Note that oj c is used instead of the actual 
circular velocity v c = Ru c since the errors of oo c and R are not correlated, while those of v c 
and R are. The total acceleration is then given by d(p t ot/dR = 

For the baryonic component, we implement a wide range of alternative observation- 


based distributions for the stellar bulge (5 

itanek et al. 

1997 

Zhao 

1996 

Bissantz & Gerhard 

2002 Lopez-Corredoira et al. 2007| Vanhollcbeke et al. 1 

1009; Robin et al. 2012), stellar 

disc(s) (Han & Gould 

2003; Calchi Novati & Mancini 

2011 

de Jong et al. 2010 Juric et al. 

2008 Bovy & Rix 20L‘ 

1) and gas (Ferriere 

1998; Moskalenko et al. 

2002 

). The bulge models 


comprise different parameterizations for the Galactic bar and are normalised to microlensing 


optical depth data (Popowski et al. 2005) using a procedure thoroughly described in Iocco 


et al. (2011). The discs bracket a variety of profiles fitted to photometric observations and 


are calibrated to the latest measurement of the local total stellar surface density (Bovy & 


Rix 2013). The two alternative gas models adopted (Ferriere 1998 Moskalenko et al. 2002) 


consist of molecular, atomic and ionised phases whose spatial distributions have been traced 
by wide surveys of different spectral lines (mainly 21 cm and CO); we have checked that using 


the more massive HI disc of Kalberla & Dedes (2008) has no significant impact on our results. 


The contribution of each component to the rotation curve is computed through multipole 
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expansion (Binney & Tremaine 2008) and the statistical uncertainty on its normalisation is 


propagated accordingly. By summing in quadrature the contribution of the three components 
(bulge, disc, gas) in all possible combinations, we obtain a compilation that covers the entire 
range of morphologies available in the literature, thus bracketing the baryonic contribution 
to the total acceleration, dcpb/dR = co%R. It is crucial to notice here that, to the best of 
our knowledge, our library of observation-inferred baryonic distributions includes virtually 
all morphologies available in the literature, which allows us to adopt the resulting spread 
as an estimate of the systematic uncertainty on baryonic modelling rather than use a single 
model with unknown systematics. 

The core aim of this work is to subtract C0b from co c . Defining dcpdm/^R = w % m R > Hie 
decomposition of the potential implies 


co 


2 

dm 



(i) 


under the assumption that the discrepancy between the observed rotation curve and that 
expected from the distribution of observed baryons is caused by an underlying dark matter 
component. The inferred residuals uo 2 drn and the corresponding uncertainties (propagated 
from both oo c and tOb) are shown in Fig. [l] for a baryonic model comprising a specific bulge 
(Stanek et al. 1997[ ), disc (Bovy & Rix 2013[ ) and gas ( |Ferriere 1998), chosen for representative 
purposes as it lies close to the median value of the baryonic envelope at all R (see next 
section). All objects within 2.5kpc from the Galactic centre are omitted to avoid tracers 
with non-circular orbits. The residuals are consistently above zero and grow towards the 
centre. These data trace dcftdm/dR and not directly the dark matter density distribution 
Pdm■ However, the radial slope of co dm = v dm /R 2 does contain information about pd m ■ For 
simplicity, let us take a spherically symmetric dark matter component. Then, one obtains 
the well-known relation v dm = GMd m (<R) / R, which can be easily solved for the density in 
a spherical shell of radius R: 


Pdm — 


47 tG 


3 co, 


dm 


+ r ^r 


UJ 


dm 


4v tG 


3 + 


din col 


dm 


din A 


( 2 ) 


where the last step is strictly valid only for co dm > 0. The same result follows directly from 
the Poisson equation for a spherical potential. Effectively, any deviation from the scaling 
co dm oc R~ 3 indicates the presence of dark matter and the magnitude of such deviation is a 
measure of its density at radius R. For the emblematic case of a flat rotation curve Vd m — 
v c = const, the usual scaling p dm oc R~ 2 is recovered. Eq. ([2]) is our master formula to extract 
the (spherical) dark matter profile directly from the data. Notice that (i) no assumption has 
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been made about the functional form of p dm (R), and (ii) in principle it is possible to find 
the equivalent of Eq. Q for non-spherical geometries by solving d<j) dm {p<im) /d-R = u dm R. 



Fig. 1.— The dark matter contribution to the rotation curve of our Galaxy. The blue data 
points show the dark matter residuals co dm as inferred from the latest compilation of rotation 


curve measurements and for the representative baryonic model (Stanek et al. 1997 Bovy 


& Rix 2013; Ferriere 1998). The red data points display the residuals after applying the 


binning procedure described in the text. The horizontal bars in the binned residuals are 
solely to indicate the length of the radial bins. Note that the dark matter contribution to 
the actual rotation curve v c reads Vd m = Ru d m- 


The determination of the profile requires an estimate of c o dm and its radial slope. The 
first step adopted is to bin the data, which comprises 2687 individual measurements in the 
range R = 2.5 — 25kpc. We start by setting up 18 linearly spaced intervals in this range 
and then merge adjacent intervals as necessary to have at least five measurements per bin 
and a mean uncertainty in R smaller than the bin half-width. Next, we compute the simple 
average of u dm in each bin and take for its uncertainty the quadrature of the mean of u dm 
uncertainties and the standard deviation of the central values. Finally, data points more 
than five sigma away from the u> dm average are excluded and the procedure is repeated 
until convergence. Typically, about 12 outliers are excluded from the initial set, and the 
remaining measurements are distributed into 10 radial bins. The number of measurements 
per bin range from 19 (for the outermost bin) to around 575. The resulting binned residuals 
are shown in red in Fig. [l] and are used to set c u dm in Eq. (J2|. 
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The next step is to determine the slope of cu^ m , namely the second term in Eq. (J2|. 
Instead of using the average residuals explained above (which would lead to unnecessary 
correlations and overestimated uncertainties), we use the individual data points within each 
bin. Specifically, in what we call our default method and for which the final results are shown, 
a weighted straight-line fit of the points in each bin is performed to estimate duj^ m / d R and 
the corresponding uncertainty. In order to validate the radial slope estimates of our default 
method, we have implemented three alternative procedures. In the first, a proxy of the 
default method, the measurements within each bin were fitted with a power law (rather than 
a straight line) to determine the logarithmic slope in the second line of Eq. ([2]). Separately, 
the slope duj^ m /dR has been estimated as a differential between adjacent bins with two 
different methods. The first relies simply on the difference between the average values of 
u>d m across adjacent bins. The second is slightly more involved: dui^ m /dR has been computed 
for each single data point in every bin as an average of the differentials obtained with all 
the data points of the previous bin. The slope doj^ m /dR assigned to the bin is then the 
average of all single-point values, and uncertainties are computed as a spread around that 
central value. The last two methods (based on the adjacent bin differential) present the 
disadvantage of having highly correlated uncertainties. However, the central values for the 
slope obtained with all four methods are in remarkable agreement, as discussed in the next 
section. 


3. Results 

We now have all the necessary ingredients to determine the dark matter profile using 
Eq. ([2]). This was done across Galactocentric radii R = 2.5 —25kpc for each baryonic model. 
We thus obtain an envelope for pd m which encompasses both the statistical uncertainty arising 
from the residual co^ m and its slope, and the systematic uncertainty due to our ignorance of 
the actual morphology of stars and gas in the Galaxy. Fig. [2] shows the determination of 
the profile obtained by applying the default method to compute the radial slope du^ m /dR. 
The error bars represent the one sigma uncertainties for the representative baryonic model; 
note that these uncertainties result from the propagation of errors on the rotation curve u c , 
the normalisation of the baryonic component ojb, the Galactocentric radius R and the slope 
d u>dm / dR. The grey region encompasses the one sigma determination for all baryonic models 
implemented. 

Before discussing Fig. [2j a few comments are in order regarding the robustness of our 
findings. Firstly, we note that the four methods devised to estimate the radial slope (see 
previous section) are all compatible with each other at the one sigma level, thus speaking for 
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Fig. 2.— The spherical dark matter profile of our Galaxy as inferred directly from observa¬ 
tions. The red data points represent the one sigma measurement of the profile for the radial 
binning shown in Fig. [I] and the representative baryonic model (Stanek et al. 1997; Bovy & 


Rix 2013 Ferriere 1998). The grey region encompasses the results for all baryonic models, 


including the corresponding one sigma uncertainties. Overplotted in black are commonly 
used dark matter profiles, namely Navarro-Frenk-White with scale radius r s = 20 kpc (solid 
line), Einasto with r s = 20kpc and slope parameter a = 0.17 (dashed) and cored isothermal 
with r s = 5 kpc (dotted), all normalised to a local dark matter density of 0.4GeV/cm 3 . For 
reference, 0.38GeV/cm 3 = 0.01M o /pc 3 . 


the consistency of the presented profile determination. Moreover, we have explicitly checked 
that the results are solid against the choice of the radial binning and the non-exclusion 
of outlier data points. Notice that throughout this work R 0 = 8 kpc and v 0 = 230km/s; 
adopting different values does not qualitatively change our conclusions. 

The results in Fig. [2] present several remarkable features. In order to address these 
properly, it is important to point out that the Erst piece in Eq. (|2]) dominates over the second 
one across the whole range of Galactocentric radii addressed here. The slope of the dark 
matter residuals is therefore sub-leading (but not negligible) in the determination of pd m - We 
first comment on the magnitude of the uncertainties shown in Fig. [2] In the innermost bins, 
both the uncertainty for a single baryonic model and the dispersion due to baryonic modelling 
are large. This is because baryons dominate the gravitational potential below about 5 kpc. 






























On the one hand, the little room left for a dark matter contribution to the rotation curve 
prompts the considerable uncertainty of the reconstructed pd m for each baryonic model. On 
the other hand, the leading role of baryons in this region gives weight to a broad dispersion in 
the contribution of different morphologies of the inner Galaxy, which in turn shows up in the 
extended size of the grey region below about 5 kpc. Such effect is mitigated between 6 and 
lOkpc, where the baryons give a decreasingly important contribution to the gravitational 
potential. In that intermediate region, the mild uncertainties reported are dominated by 
the dispersion of the rotation curve data. This dispersion then grows towards larger radii, 
causing the fluctuations seen above lOkpc. 


The uncertainties in the innermost regions will eventually be improved with data soon to 
be provided by Gaia (de Bruijne 2012), whose dramatic impact on the census of the Galaxy 
will very likely help reduce the spread on the current models of the inner Milky Way. For 
larger Galactocentric radii, an increase in number and precision of rotation curve measure¬ 
ments (or the use of alternative kinematic tracers) would also improve the reconstruction 
of the dark matter profile. Interestingly enough, although our method is not optimised to 
measure the dark matter density in the solar neighbourhood, we do find a density at R ~ i? 0 
close to the usual values 0.3 — 0.5GeV/cm 3 obtained by both global (Dchnen & Binney 


1998 

Sofue et ah 2C 

09 Catena & Ullio 2010 

Weber & de Boer 2010 Iocco et ah 2011 

Nesti & Salucci 

2013 

Bovy & Rix 

2013 

) and local ( 

Salucci et ah 

2010 

Smith et ah 

2012 


Garbari et ah 2012 Zhang et al. 2013 Read 2014) methods. Notice in particular that our 


results are compatible with the findings of Salucci et ah (2010) even though we use no con¬ 
straints on the Oort’s constants to fix the local slope of the rotation curve. Furthermore, as 
shown in Fig. [2j our reconstructed profile is in agreement with those inferred from numerical 
simulations (Navarro et ah 1996 Merritt et ah 2006), but current uncertainties hinder any 
discrimination power between different radial behaviours. In principle, it would be possible 
to shrink the reported uncertainties at the expense of forcing a generic functional form for 
the dark matter profile (e.g., a monotonicity prior), but we refrain to do so here in order 
not to spoil the innovative feature of our technique. Our minimal approach also allows for 
future observational tests of theoretical priors, thus making it a powerful diagnostic tool. 


4. Conclusion 

It is truly remarkable that, despite decades of observations and theoretical progress, the 
distribution of dark matter in the Milky Way remains largely unconstrained. The situation 
is particularly problematic towards the inner Galaxy, where baryons dominate the gravita¬ 
tional potential and for which any solid improvement would have a remarkable impact on 


















































































astroparticle physics and cosmology. In this context, we have for the first time implemented 
a method to reconstruct the Galactic dark matter profile directly from observations. The 
method requires no assumption on the form or shape of the profile, unlike all previous tech¬ 
niques applied to the Milky Way. Our findings - obtained using the most recent kinematic 
data and baryonic models - are in good agreement with numerical simulations and with 
both local and global measurements of the local dark matter density. These results can be 
improved both on the observational side (e.g., by including future Gaia data) and on the 
theoretical side (e.g., by relaxing the assumption of spherical symmetry). The proposed 
technique, complementary and competitive to others in the literature, represents a step for¬ 
ward towards achieving a more accurate description of the dark matter distribution in our 
Galaxy. 
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